Il s’agira de mener une série de tests statistiques sur l’archive du même intitulé. Cet archive (au format CSV) contient des entrées relatives aux faits divers causés par des armes à feu qui ont eu lieu aux Etats Unis entre le 1er janvier 2013 et le 31 mars 2018. Il s’agira de montrer l’évolution d’un certain nombre de variables pendant la période concernée, d’en calculer la moyenne et l’écart type (quand cela est possible). L’ensemble des variables qui nous intéressent contient : — nombre de blessés ; — nombre de morts ; — nombre de malfaiteurs ; — age des malfaiteurs ; — état dans lequel le fait divers a eu lieu.
Le dataset ne contient pas la variable nombre de malfaiteurs, nous ne pouvons donc pas la traiter.
Le sujet nous demande ici d’étudier une base de donnée regroupant diverses informations sur des faits causés par les armes aux Etats Unis entre le 1er Janvier 2013 et le 31 Mars 2018. Le but de cet étude est d’observer l’évolution de certaines variables au fur et a mesure du temps. C’est variables concerne le nombre de blessés, le nombre de morts, le nombre de malfaiteurs, leur age et l’état dans lequel le fait divers a eu lieu.
Les données de cette étude ont été récoltés par une société à but non lucratif créée en 2013 pour fournir un accès au public des informations précises sur la violence liée aux armes à feu aux États-Unis. Cette société a collecté et vérifié l’exactitude des informations. La base de donnée comporte pas moins de 260 000 incidents de violence armée survenue entre entre le 1er Janvier 2013 et le 31 Mars 2018 aux Etats Unis. Les données sont triées par date (ordre croissant). Chaque donnée comporte une variété d’informations, dont un id, une date, l’état ou s’est déroulé le crime, la ville, l’adresse, le nombre de morts, le nombre de blessés, le type d’arme utilisé, le nombre d’armes a feu et même l’age et le genre des participants.
library(knitr)
library(dplyr)
library(readr)
library(ggplot2)
library(tibble)
library(stringr)
library(gridExtra)
library(scales)
library(lubridate)
library(tibble)
library(purrr)
library(splitstackshape)
library(PerformanceAnalytics)
library(tidyr)
library(corrplot)
library(lubridate)
On commence par créer une dataframe (df) en indiquant que tout les NA et les espaces sont des NA :
df = read.csv('/Users/ALEXIS/Downloads/gun-violence-data_01-2013_03-2018.csv', header=TRUE,stringsAsFactors = FALSE, na.strings=c("NA", ""))
Nous supprimons les colonnes inutiles dans une nouvelle dataframe (df1):
df1 <- subset(df, select = -c(incident_url_fields_missing,location_description,notes,incident_url,incident_characteristics,longitude,latitude,sources, source_url))
#Valeurs manquantes
apply(df1, 2, function(col)sum(is.na(col))/length(col))
## incident_id date state
## 0.00000000 0.00000000 0.00000000
## city_or_county address n_killed
## 0.00000000 0.06883013 0.00000000
## n_injured congressional_district gun_stolen
## 0.00000000 0.04983373 0.41513370
## gun_type n_guns_involved participant_age
## 0.41493760 0.41493760 0.38509327
## participant_age_group participant_gender participant_name
## 0.17573234 0.15171251 0.51007397
## participant_relationship participant_status participant_type
## 0.93418643 0.11526346 0.10373544
## state_house_district state_senate_district
## 0.16176771 0.13491073
Nous créons une base de données ne comportant que les variables qui nous intéressent. Cette base de donnée contient beaucoup de valeurs manquantes. Alors nous mesurons en pourcentage la qualité des données de chaque colonnes pour ensuite supprimer celles qui ont plus de 30% de valeurs manquantes :
#We'll delete all the columns with more than 40% of missing values
df1 <- subset(df1, select = -c(gun_stolen,gun_type,n_guns_involved,participant_name,participant_relationship))
apply(df1, 2, function(col)sum(is.na(col))/length(col))
## incident_id date state
## 0.00000000 0.00000000 0.00000000
## city_or_county address n_killed
## 0.00000000 0.06883013 0.00000000
## n_injured congressional_district participant_age
## 0.00000000 0.04983373 0.38509327
## participant_age_group participant_gender participant_status
## 0.17573234 0.15171251 0.11526346
## participant_type state_house_district state_senate_district
## 0.10373544 0.16176771 0.13491073
Fondamentalement, nous devrions supprimer participant_age, car cette colonne contient 38,5 % de données manquantes, Mais nous devons le garder pour l’étude.
Nous effaçons toutes les valeurs manquantes, nous pourrions les remplacer par la méthode Miss FAMD, mais comme nous ne faisons pas de machine learning ce n’est donc pas nécessaire.
df1=na.omit(df1)
Nous formatons certaines données (celles comportant des : ||) afin d’avoir des valeurs numériques et une meilleure lisibilité:
#We'll change some data for the readibility of them
age=cSplit(df1,c("participant_age"),sep="||",direction="wide",drop=FALSE)
age$age=gsub(".*::","",age$participant_age)
age$age=as.numeric(age$age)
df1$participant_age=age$age
type=cSplit(df1,"participant_status",sep="||" ,direction="wide",drop=FALSE)
type$participant_status=gsub(".*::","",type$participant_status)
df1$participant_status=type$participant_status
type=cSplit(df1,"participant_age_group",sep="||",direction="wide",drop=FALSE)
type$participant_age_group=gsub(".*::","",type$participant_age_group)
df1$participant_age_group=type$participant_age_group
type=cSplit(df1,c("participant_gender"),sep="||",direction="wide",drop=FALSE)
type$participant_gender=gsub(".*::","",type$participant_gender)
df1$participant_gender=type$participant_gender
type=cSplit(df1,c("participant_type"),sep="||",direction="wide",drop=FALSE)
type$participant_type=gsub(".*::","",type$participant_type)
df1$participant_type=type$participant_type
Nous créons un nouveau dataframe(df2) à partir des colonnes de la précédente dataframe(df1) en gardant seulement les colonnes désirés pour l’étude :
#Now, after some preprocessing and the reorganization of the columns, we will kept only the
#columns needed.
#State, n_killed, n_injured, participant_ages
df2 <- subset(df1, select = c(n_killed, n_injured, state, participant_age))
head(df2)
## n_killed n_injured state participant_age
## 2 1 3 California 20
## 3 1 3 Ohio 33
## 4 4 0 Colorado 33
## 5 2 2 North Carolina 47
## 6 4 0 Oklahoma 55
## 7 5 0 New Mexico 15
Nous supprimons les NA qui restent et les valeurs peu probables(valeurs abérrantes) :
#We delete the last NA
df2=na.omit(df2)
summary(df2$participant_age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 21.00 26.00 29.67 35.00 209.00
#We can see a crazy value, we will delete it
df2 <- df2[-which(df2$participant_age==209),]
dfnum <- subset(df2, select = c(n_killed, n_injured, participant_age))
Afin de voir s’il y a des corrélations nous avons utilisé plusieurs méthodes, la première, est visuelle, à l’aide de plot, afin de vérifier s’il y a des tendances ou non. La seconde est statistique à l’aide des matrices de corrélations, nous devons préalablements formatés les données afin d’utiliser celles qui sont numériques
Ici nous produisons des graphiques, par rapport aux états et nombres de morts et aux états et nombres de bléssés
bystateinj <- df2 %>% group_by(state) %>% summarize(n = n(), ninj = sum(n_injured)) %>% tidyr::gather(key = "type", value = "num", ninj)
ggplot(bystateinj, aes(x = state, y = num, fill = type)) + geom_boxplot() + theme_bw()
bystatekill <- df2 %>% group_by(state) %>% summarize(n = n(), nkill = sum(n_killed))%>% tidyr::gather(key = "type", value = "num", nkill)
ggplot(bystatekill, aes(x = state, y = num, fill = type)) + geom_boxplot() + theme_bw()
bystate <- df2 %>% group_by(state) %>% summarize(n = n(), nkill = sum(n_killed), ninj = sum(n_injured)) %>% tidyr::gather(key = "type", value = "num", nkill, ninj)
ggplot(bystate, aes(x = state, y = num, fill = type)) + geom_boxplot() + theme_bw()
Ici nous faisons une matrice de corrélation de Pearson, qui nous montre l’absence de corrélation
summary(df2$participant_age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 21.00 26.00 29.67 35.00 101.00
cor(dfnum, method = c("pearson"))
## n_killed n_injured participant_age
## n_killed 1.0000000 -0.1601249 0.1224553
## n_injured -0.1601249 1.0000000 -0.1073301
## participant_age 0.1224553 -0.1073301 1.0000000
mcor <- cor(dfnum)
mcor
## n_killed n_injured participant_age
## n_killed 1.0000000 -0.1601249 0.1224553
## n_injured -0.1601249 1.0000000 -0.1073301
## participant_age 0.1224553 -0.1073301 1.0000000
col<- colorRampPalette(c("blue", "white", "red"))(20)
heatmap(x = mcor, col = col, symm = TRUE)
Ici nous faisons un test de corrélation
respear1nkni<-cor.test(dfnum$n_killed,dfnum$n_injured, method="pearson")
respear1nkni
##
## Pearson's product-moment correlation
##
## data: dfnum$n_killed and dfnum$n_injured
## t = -54.046, df = 111000, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1658514 -0.1543875
## sample estimates:
## cor
## -0.1601249
respear1nkage<-cor.test(dfnum$n_killed,dfnum$participant_age, method="pearson")
respear1nkage
##
## Pearson's product-moment correlation
##
## data: dfnum$n_killed and dfnum$participant_age
## t = 41.107, df = 111000, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1166566 0.1282457
## sample estimates:
## cor
## 0.1224553
Comme nous pouvons le voir il n’y a aucun correlation entre les subsets
Encore une fois nous pouvons remarquer l’absence de corrélations fortes.
Nous commençons par séléctionner la partie de la data que nous avons besoin (la période entre le 1er janvier 2013 et le 31 décembre 2017 et celle entre le 1 janvier 2018 et le 31 mars 2018) :
#Prenons les datas df1 avec les colonnes utiles
data <- subset(df1, select = c(n_killed, n_injured, state, participant_age, date))
#et prenons les datas recueillis en 2018 (01-01-2018 -> 31-12-2018)
data_2018 <- subset(data, data$date >= as.Date('2018-01-01') & data$date <= as.Date('2018-03-31'))
#Pour la période suivante : 01-01-2013 -> 31-12-2017
data_With_Periode <- subset(data, data$date >= as.Date('2013-01-01') & data$date <= as.Date('2017-12-31'))
#Etudions plus en détail ce jeu de données
summary(data_With_Periode)
## n_killed n_injured state participant_age
## Min. : 0.0000 Min. : 0.0000 Length:106480 Min. : 0.00
## 1st Qu.: 0.0000 1st Qu.: 0.0000 Class :character 1st Qu.: 21.00
## Median : 0.0000 Median : 0.0000 Mode :character Median : 26.00
## Mean : 0.3777 Mean : 0.5534 Mean : 29.63
## 3rd Qu.: 1.0000 3rd Qu.: 1.0000 3rd Qu.: 35.00
## Max. :50.0000 Max. :53.0000 Max. :209.00
## NA's :1597
## date
## Length:106480
## Class :character
## Mode :character
##
##
##
##
#a voir si on peut améliorer cette partie
#Considérons l'échantillon suivant : data_2018
summary(data_2018)
## n_killed n_injured state participant_age
## Min. : 0.0000 Min. : 0.0000 Length:6120 Min. : 0.00
## 1st Qu.: 0.0000 1st Qu.: 0.0000 Class :character 1st Qu.:21.00
## Median : 0.0000 Median : 0.0000 Mode :character Median :27.00
## Mean : 0.3956 Mean : 0.4667 Mean :30.41
## 3rd Qu.: 1.0000 3rd Qu.: 1.0000 3rd Qu.:37.00
## Max. :17.0000 Max. :17.0000 Max. :94.00
## date
## Length:6120
## Class :character
## Mode :character
##
##
##
Si nous comparons les résultats, nous remarquons les choses suivantes :
H0 NB MORT : Le nombre de morts obtenus en 2018 a permis de modifier la moyenne de mort HO NB BLESSE : Le nombre de blesses obtenus en 2018 a permis de modifier la moyenne de blesses HO AGE PARTICIPANT : L’age des participants releves en 2018 a permis de modifier la moyenne d’age
Pour les variables relatives aux nombres de mort, blesses et l’age des protagonistes : nous effectuerons un test bilatérale afin de valider ou non H0
Nous etablissons nos fonctions au préalable ainsi elles faciliterons le traitement et ce sera plus lisible pour la suite.
test <- function(mu, n, xbar, sigma) {
res = abs(xbar-mu)/(sigma/sqrt(n))
return(res)
}
p_value <- function(alpha, ddl) {
pval = qt(alpha, df=ddl, lower.tail = FALSE)
return(pval)
}
sigma_value <- function(data) {
res = var(data)
return(res)
}
#taille de l'echantillon
n <- length(data_2018)
# On trouve notre valeur de comparaison (P-value) avec alpha = 5%
tPVAL <- p_value(0.05, n-1)
#d'après les resultats obtenus, on a :
mu_NbMort <- mean(data_2018$n_killed)
Xbar_NbMort <- mean(data_With_Periode$n_killed)
sigma_NbMort <- sigma_value(data_2018$n_killed)
# On applique notre test avec les données obtenues
t_NbMort <- test(mu_NbMort, n, Xbar_NbMort, sigma_NbMort)
# On compare et on regarde le résultat
if(t_NbMort >= -tPVAL && t_NbMort <= tPVAL) {
cat("pour alpha = 5%, on a t_nbMort (", t_NbMort ,") appartenant à l'intervalle [",
-tPVAL , ";", tPVAL, "]\n-> Ainsi, on ne peut pas refuser H0\n")
} else {
cat("pour alpha = 5%, on a t_nbMort (", t_NbMort ,") n'appartenant pas à l'intervalle [",
-tPVAL , ";", tPVAL, "]\n-> Ainsi,on rejette H0\n")
}
## pour alpha = 5%, on a t_nbMort ( 0.09824325 ) appartenant <U+00E0> l'intervalle [ -2.131847 ; 2.131847 ]
## -> Ainsi, on ne peut pas refuser H0
#d'après les resultats obtenus, on a :
mu_NbBlesses <- mean(data_2018$n_injured)
Xbar_NbBlesses <- mean(data_With_Periode$n_injured)
sigma_NbBlesses <- sigma_value(data_2018$n_injured)
# On applique notre test avec les données obtenues
t_NbBlesses <- test(mu_NbBlesses, n, Xbar_NbBlesses, sigma_NbBlesses)
# On compare et on regarde le résultat
if(t_NbBlesses >= -tPVAL && t_NbBlesses <= tPVAL) {
cat("pour alpha = 5%, on a t_nbMort (", t_NbBlesses ,") appartenant à l'intervalle [",
-tPVAL , ";", tPVAL, "]\n-> Ainsi, on ne peut pas refuser H0\n")
} else {
cat("pour alpha = 5%, on a t_nbMort (", t_NbBlesses ,") n'appartenant pas à l'intervalle [",
-tPVAL , ";", tPVAL, "]\n-> Ainsi,on rejette H0\n")
}
## pour alpha = 5%, on a t_nbMort ( 0.3726553 ) appartenant <U+00E0> l'intervalle [ -2.131847 ; 2.131847 ]
## -> Ainsi, on ne peut pas refuser H0
#d'après les resultats obtenus, on a :
mu_Age <- mean(data_2018$participant_age)
Xbar_Age <- 29.63 # mean(data_With_Periode$participant_age) la cmd ne marche pas ... à voir
#tmp : on la saisie à la main
sigma_Age <- sigma_value(data_2018$participant_age)
# On applique notre test avec les données obtenues
t_Age <- test(mu_Age, n, Xbar_Age, sigma_Age)
# On compare et on regarde le résultat
if(t_Age >= -tPVAL && t_Age <= tPVAL) {
cat("pour alpha = 5%, on a t_nbMort (", t_Age ,") appartenant à l'intervalle [",
-tPVAL , ";", tPVAL, "]\n-> Ainsi, on ne peut pas refuser H0\n")
} else {
cat("pour alpha = 5%, on a t_nbMort (", t_Age ,") n'appartenant pas à l'intervalle [",
-tPVAL , ";", tPVAL, "]\n-> Ainsi,on rejette H0\n")
}
## pour alpha = 5%, on a t_nbMort ( 0.01011114 ) appartenant <U+00E0> l'intervalle [ -2.131847 ; 2.131847 ]
## -> Ainsi, on ne peut pas refuser H0
Si l’on prend en compte aussi le mois de l’année dans lequel le fait divers a été commis, est-ce qu’il y a une corrélation forte entre le nombre de fait divers et le mois de l’année ? Quels conclusions en tirez-vous ?
Ici, nous devons créer deux nouvelles colonnes afin de traiter le mois et l’année. Mais avant ça, il faut faire du preprocessing sur les données afin d’harmoniser la colonne date et de la rendre utilisable. On va donc la formater sous forme année, mois, jour.
df3 <- subset(df1, select = c(n_killed, n_injured, state, participant_age, date))
df3$date <- ymd(df3$date)
str(df3$date)
## Date[1:112600], format: "2013-01-01" "2013-01-01" "2013-01-05" "2013-01-07" "2013-01-07" ...
A la suite de ça, nous mappons les variables de mois entre 1 et 12 afin de réutiliser des matrices de corrélations On supprime les dernieres valeurs “NA” qui se sont rajoutés.
df3=na.omit(df3)
df3$month <- month(df3$date, label=TRUE)
df3$year <- year(df3$date)
df3$month <- factor(df3$month, levels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"), labels = c(1,2,3,4,5,6,7,8,9,10,11,12))
Produisons des graphiques afin de vérifier la corrélation entre les variables. Dans le dernier graphique affiché, on peut voir qu’il n’ya pas de tendance entre les valeurs, la seule tendance que l’on peut remarquer est que le mois de juillet est le plus meurtrier. Certainement avec le 4 Juillet
plotly::ggplotly(df3 %>% select(month) %>% count(month) %>%
ggplot(aes(x=month, y=n)) + geom_bar(stat='identity', fill='red') +
scale_y_continuous(labels=comma) +
labs(x='Months', y='Number of incidents', title='Incidents by Month'))
plotly::ggplotly(df3 %>% select(year, month) %>% group_by(year) %>% count(month) %>%
ggplot(aes(x=month, y=n)) + geom_bar(stat='identity', fill='red') +
scale_y_continuous(labels=comma) + facet_grid(.~year) +
labs(x='Months', y='Number of incidents', title='Incidents by Month by year'))
Pour finir on utilise une nouvelle fois des matrices de corrélations afin de vérifier si corrélation il y a ou non. On doit donc transformer les variables de type string en variable numérique Et comme vous pouvez le voir, il n’y en a pas.
df3num <- subset(df3, select = c(n_killed, n_injured, participant_age, month, year))
class(df3$month)
## [1] "ordered" "factor"
df3num$month=as.numeric(df3num$month)
df3num$year=as.numeric(df3num$year)
df3num<-na.omit(df3num)
cor(df3num, method = c("pearson"))
## n_killed n_injured participant_age month
## n_killed 1.000000000 -0.16013017 0.1224744858 -0.0015462327
## n_injured -0.160130170 1.00000000 -0.1073241728 0.0147094425
## participant_age 0.122474486 -0.10732417 1.0000000000 -0.0008562048
## month -0.001546233 0.01470944 -0.0008562048 1.0000000000
## year -0.009938460 -0.01986516 0.0105802256 -0.1955934250
## year
## n_killed -0.00993846
## n_injured -0.01986516
## participant_age 0.01058023
## month -0.19559343
## year 1.00000000
mcor1 <- cor(df3num)
mcor1
## n_killed n_injured participant_age month
## n_killed 1.000000000 -0.16013017 0.1224744858 -0.0015462327
## n_injured -0.160130170 1.00000000 -0.1073241728 0.0147094425
## participant_age 0.122474486 -0.10732417 1.0000000000 -0.0008562048
## month -0.001546233 0.01470944 -0.0008562048 1.0000000000
## year -0.009938460 -0.01986516 0.0105802256 -0.1955934250
## year
## n_killed -0.00993846
## n_injured -0.01986516
## participant_age 0.01058023
## month -0.19559343
## year 1.00000000
col<- colorRampPalette(c("blue", "white", "red"))
corrplot(mcor1, type="upper", order="hclust", tl.col="black", tl.srt=45)
chart.Correlation(df3num, histogram=TRUE, pch=19)
Ici nous allons refaire des tests de corrélations afin de vérifier si corrélation il y a.
respear3nkni<-cor.test(df3num$n_killed,df3num$n_injured, method="pearson")
respear3nkni
##
## Pearson's product-moment correlation
##
## data: df3num$n_killed and df3num$n_injured
## t = -54.048, df = 111000, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1658567 -0.1543928
## sample estimates:
## cor
## -0.1601302
respear3nkage<-cor.test(df3num$n_killed,df3num$participant_age, method="pearson")
respear3nkage
##
## Pearson's product-moment correlation
##
## data: df3num$n_killed and df3num$participant_age
## t = 41.114, df = 111000, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1166758 0.1282648
## sample estimates:
## cor
## 0.1224745
respear3niage<-cor.test(df3num$n_injured,df3num$participant_age, method="pearson")
respear3niage
##
## Pearson's product-moment correlation
##
## data: df3num$n_injured and df3num$participant_age
## t = -35.965, df = 111000, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1131355 -0.1015055
## sample estimates:
## cor
## -0.1073242
respear3nkyr<-cor.test(df3num$n_killed,df3num$year, method="pearson")
respear3nkyr
##
## Pearson's product-moment correlation
##
## data: df3num$n_killed and df3num$year
## t = -3.3113, df = 111000, p-value = 0.0009288
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.015820301 -0.004055931
## sample estimates:
## cor
## -0.00993846
respear3nimth<-cor.test(df3num$n_injured,df3num$month, method="pearson")
respear3nimth
##
## Pearson's product-moment correlation
##
## data: df3num$n_injured and df3num$month
## t = 4.9012, df = 111000, p-value = 9.536e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.00882744 0.02059043
## sample estimates:
## cor
## 0.01470944
C’est encore loin de 1 et de -1 il n’y a donc pas de corrélations.